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Forty years ago Richardson showed that the eigenstates of the pairing Hamihonian with constant 
interaction strength can be calculated by solving a set of non-linear coupled equations. However, 
in the case of Fermions these equations lead to singularities which made them very hard to solve. 
This letter explains how these singularities can be avoided through a change of variables making 
the Fermionic pairing problem numerically solvable for arbitrary single particle energies and degen- 
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Introduction 



Exactly solvable models serve as a guideline for un- 
derstanding the properties of correlated many-body sys- 
tems. Already in the sixties, Richardson solved the eigen- 
problem of a constant pairing interaction in a set non- 
degenerate single-particle levels for Fermion Q and Bo- 
son '5| systems. However, it turned out that in the 
Fermion case the solutions exhibit singularities which 
are hard to treat numerically 0, Q. Recently, exactly 
solvable pairing models have gained new attention [3, 
with applications to nanometalic grains [6j (for a re- 
view see 13), Bose Einstein Condensates ^ and nuclear 
physics It was first shown that the pairing model 
was integrable by finding the complete set of commuting 
integrals of motions |lO|, and subsequently, three new 
families of fully integrable and exactly solvable models 
[lll | giving rise to a large class of pairing Hamiltonians 
with non-uniform matrix elements jl2l | were presented. 
These models are exactly solvable, except for the sin- 
gularities occurring in Fermion systems for some critical 
values of the pairing strength. T his p roblem, in spite of 
some early attempts to cure it 0, 0| , precluded for over 
forty years the use of these exactly solvable models for a 
wide range of applications, ranging from condensed mat- 
ter to nuclear physics. Moreover, the recent developed 
extensions of the exact solution to pairing Hamiltonians 
including the isospin degree of freedom (0(5) pairing) 
[Tsl with promising applications io N Z nuclei and 
high Tc superconductivity, suffer from the same kind of 
singularities. 

In this letter show how the Richardson equations can 
be solved numerically, avoiding the singularities, through 
an appropriate change of variables. This procedure pro- 
vides a fast an accurate way to solve the equations for a 
constant pairing interaction. The method can be applied 
as well to more general exactly solvable Hamiltonians as- 
sociated with a coupled set of non- linear equations of the 
Richardson type. 

The exactly solvable pairing Hamiltonian has the fol- 



lowing form: 
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with Cj any set of singic-particle energies and g the pair- 
ing interaction strength. The iV-pair eigenstates of this 
Hamiltonian have the form 
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where |0) is a state without paired particles (a Racah 
quasispin vacuum state, fl^). The corresponding energy 
E{{x}) is given by 



E{{x}) = {0\H\0)+Y,- 



(3) 



The complex variables be found by solving a set 

of non- linear equations: 
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for a — 1, . . . , N. The parameters dj depend on the level 
degeneracies and on the structure of the vacuum state 
|0). For Fermions, dj = '^^ , where Q.j is the pair de- 
generacy (for the nuclear shell model VLj = j + 1/2 ), and 
i/j is the seniority of the level j. Note that dj < be- 
cause of the Pauli principle. Solving this set of nonlinear 
equations solves the eigenproblem for the pairing Hamil- 
tonian Eq.QJ. Unfortunately, the algebraic solution be- 
comes numerically unstable at certain critical values of 
the interaction strength [^Q . This is caused by singular- 
ities in the first and second terms in Eq.Q), when some 
of the variables Xa are approaching the value 2ej . One 
can understand this from the electro mag netic analogy for 
the exactly solvable pairing model |l7l |: in the Fermion 
case, the single-particle levels Cj and the variables Xa cor- 
respond to opposite charges. Therefore a group of vari- 
ables can cluster around a single-particle level in such a 
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FIG. 1: Real and imaginary part of the variables Xa for the 
model described in the text 



way that for each of the variables the repulsive charge of 
the other variables is compensated for by the attractive 
charge of the single-particle level. These singularities oc- 
cur for Fermions in double or multiple degenerate levels. 
In the case of doubly degenerate equidistant levels these 
singularities can be handled for the ground state , but 
the problem arises again in the treatment of the excited 
states 18]. In previous calculations it was necessary to 
tune the solutions by hand as soon as they approach a 
singularity. No general solution method was known until 
now. Figure ^ illustrates the behavior of the variables 
in case of multiply degenerate levels (see below for the 
details of the model). 

A general approach to solve the Richardson equa- 
tions Eq.|01), starts from an approximate solution in the 
weak-interaction limit (see below). Then this solution 
is evolved adiabatically up to the desired interaction 
strength by gradually increasing the value of g. At each 
step in g, the previous solution has to be updated. For 
dealing with the singularities, it is useful to notice first 
that the variables Xa in Eq.Q are most sensitive to the 
other variables nearby. Therefore one can divide the set 
of variables Xa into several clusters of variables, grouped 



around different single-particle levels. By solving the 
equations for each cluster separately, one can obtain a so- 
lution for the whole system iteratively. A practical way 
to organize the clusters, is to link each variable to its 
nearest single-particle level 2ej, and to consider a cluster 
for each level that has variables around it. 

The question now is how to solve the equations for each 
cluster, particularly in the case of singularities. Let us 
consider the set of indices Ck of the Nk variables that 
cluster around a level 2efc. Then one can consider the 
equations 
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The function Fk{x) describes the influence of the other 
levels and the variables of the other clusters on the vari- 
ables in the cluster Ck ■ Because of the way the clusters 
are set up, the function Fk{x) will be a smooth function 
in the region around where the variables of the cluster 
are located. The singularities will occur in the first two 
terms of Eq.©. In the case that some of the variables 
in the cluster approach the value 2efc, the divergences in 
the first and the second term of Eq.© must cancel out. 
Multiplying Eq.jSl by 2efc — x^, and summing over the Uk 
variables Xa at the singular point, leads to the condition 



Uk = -2dk + 1, 



(7) 



with Uk the number of variables that actually converge to 
2efe. For Fermions the value {~2dk) corresponds to the 
pair degeneracy of the level efei because of the Pauli prin- 
ciple, no more Fermion pairs can occupy the level. Trying 
to put more pairs in that level results in a singularity. In 
fact, the structure of the ground state Eq.Q does not 
result in a forbidden occupation of the level. However, 
on expanding the wave function of Eq.Q in terms of 
the pair creation operators aj^a]^, one finds that the 
leading term cancels out because of the Fermionic anti- 
commutation rules. This translates into the numerical 
difficulties encountered in the solution of the equations. 
For Bosons dfc > 0, and hence Eq.([7|l shows that sin- 
gularities do not occur in the Bosonic case. This can 
also be understood from the electrostatic analogy: Boson 
pairs and single-particle levels have charges of the same 
sign. Therefore the variables try to avoid each other at 
all times, and singularities do not occur. 

The above procedure suggests a way to remove the 
singularities from the equations: multiplying Eq.(|SJ) with 
(2efc — Xq)^, for some power p, and summing over all vari- 
ables Xa in the cluster. The resulting equations become, 
for p > 1: 
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Rp = 0, (8) 
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with 

Sp = - xc^r (9) 

aeCk 

Rp = Y,{2ek-x^)PFk{xc.). (10) 

The compact form of Eq.® suggests that it might be 
advantageous to solve them for the new variables Sp in- 
stead of the original variables Xa ■ Note that given a set of 
variables Si, ... , S'at^ , one can easily construct the poly- 
nomial whose roots correspond to the values 2ek — Xa. 
Hence one can switch from one set of variables to the 
other. The problem comes with the quantities Rp: these 
are functions of the Xa, and it is not straightforward to 
express them as functions of the Sp. However, for a given 
set of variables x^ , one can easily evaluate the values Sp 
and Rp. Furthermore, one can evaluate the gradient ma- 
trix G, with Gi„i the derivative with respect to Sm of 
Eq.p forp = ; + l. 
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for l,m = 1, . . . , Nk. G^ can be evaluated accurately us- 
ing a special inversion algorithm for Vandermonde ma- 



trices PJIJ . Therefore one can solve the new set of equa- 
tions, Eq.JSI) for p = 2, . . . A''fc -I- 1, in the new variables 
Si, . . . Sn^ using a standard gradient technique such 
as the multi-dimensional Newton-Raphson method [20l| . 
However, in the case of a singularity, the gradient ma- 
trix becomes ill-conditioned: the diagonal elements of 
the gradient matrix are given by 



Gil — dk 
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for I = 1, . . . Nk . The diagonal element will vanish for 
the index Ig = 2{dk + N^) — 1. This will occur as soon as 
Nk > —2dk + 1, which matches the value for which sin- 
gularities can occur, see Eq.|(7I). In such a case the lower- 
triangular matrix G^ becomes singular. The other part 
of the gradient matrix, G^, is derived from the smooth 
function F^. A series expansion of Fk{x) in x will be 
dominated by the lowest orders. Therefore the elements 
of are very small for larger values of to. As a result, 
the value of Si^ can not be determined accurately from 
the set of equations Eq.©. One can avoid this problem 
by limiting the cluster sizes to at most the critical value 
A^fc = "-fc = — 2(ifc -|- 1, and by using gS-i as an unknown 
variable instead of S'tv^ • If more variables are found near 
to the same single-particle level, one can always divide 



the cluster into smaller, well-separated clusters, because 
at most Uk of the variables can approach the single- 
particle level closely. Knowing S-i and Si,. . . , Sn^-i, 
one can still straightforwardly construct the polynomial 
whose roots give the corresponding values Xa ■ Therefore 
one can easily switch between the two sets of variables. 
Furthermore gS^i behaves smoothly, even at a singular- 
ity. To set up an efhcient gradient method, it is useful to 
replace Eq.® for the last value, p = Nk + 1, by a similar 
equation obtained using p = 0: 



dkgS^i + gRo = 0, 



(15) 



with i?o = J2aeCk Pk{Xa)- 

For weak interaction strengths the function Fk{xa) is 
dominated by the constant term see Eq.®. In the 
weak interaction limit one can take Fk to be a constant. 
The resulting functions Rp take the simple form 



(16) 



Now the equations Eq.® and Ea. (|15|l can be solved 
straightforwardly to yield the variables Si, from which 
one can construct the polynomial that gives a unique set 
of variables Xa ■ The resulting eigenstate will depend on 
the size of the cluster for each of the single-particle lev- 
els. This establishes a one-to-one correlation between the 
eigenstates of the non-interacting system (g = 0) and 
the eigenstates of the weakly-interacting system. One 
can conclude that the Richardson equations are complete: 
their solutions generate all eigenstates, and there are no 
spurious solutions. 

One can obtain the solutions for strong interaction 
strengths by solving the weakly interacting case first, 
and then gradually increasing the interaction strength. 
At each step one can use the gradient method outlined 
above in order to update the solution to the new interac- 
tion strength. It is useful to adapt the stepsize in inter- 
action strength to the convergence of the iterative proce- 
dure by taking smaller steps in g when the convergence 
of the Newton-Raphson method for the variables Si be- 
comes slower, which typically occurs around the critical 
g values. One more ingredient is needed to avoid prob- 
lems with the singularities: when the interaction strength 
passes through a critical value, the variables Xa passing 
through a singularity can change from real to complex 
or vice versa. At the same time the variables Si will 
become very small, except for I = —1. Even using the 
new variables, the gradient method does not lead to the 
right solution when it has to pass through critical values 
of the interaction strength. One can avoid this problem 
by including an extrapolation step based upon the pre- 
vious solutions. This extrapolation has to be done in the 
variables Si , because they vary smoothly through the sin- 
gularities and they remain real all the time. Assume that 
converged solutions x'^ and x'^ were obtained for values 
g' and g" of the interaction strength. The variables x'^ 
are grouped into clusters. For each cluster, one evaluates 
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the variables S'^ and S'^. Because the variables Sp behave 
smoothly, even near a singularity, one can estimate the 
variables Sp for the new interaction strength g by linear 
extrapolation: 

, {g-g')S';-{g-g")S'p 

The resulting values of Sp can then be updated using the 
Newton-Raphson method or another gradient method. 
The extrapolation step avoids the problems with the sin- 
gularities and greatly improves the convergence of the 
method. 
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TABLE I: Woods-Saxon single-particle levels 



As an example, consider the level scheme listed in Ta- 
ble together with a constant pairing interaction. This 
model describes neutrons in ^^Fe; its ground state and 
finite-temperature properties have been studied using a 
Quantum Monte Carlo method pj| . The eigenstates can 
also be found through the solution of Richardson's equa- 
tions, Eg . lHI) , or through Lanczos diagonalization in a se- 
niority basis [221 ■ The full many-body space has a dimen- 
sion of the order of 10^*^, while the zero-seniority basis 
has dimension 14894. In Fig. [21 the lowest zero-seniority 
eigenvalues for this model are shown as a function of 
the interaction strength, calculated using the Lanczos 
method and using the method explained above (only 50 
Lanczos iterations were used). To calculate the ground 
state, a straightforward implementation of the Newton- 
Raphson method for the original equations, Eq.l^J, works 
well up to an interaction strength of ~ 0.2. It turns 
out that a singularity occurs around the 1^3/2 level at a 
value oi g = 0.245. A Fortran-95 computer program was 
written based upon the procedure outlined above. It was 
able to solve the equations for all interaction strengths 
in a matter of seconds. Figure [3| shows the behavior of 
the three x variables that cluster around the level 
as a function of g. In Fig. 0| one can see that the corre- 
sponding variables Si behave much more smoothly. More 
singularities occur around other levels at higher interac- 
tion strengths, as is shown in Fig. ^ The solution of the 
Richardson equations is faster than the Lanczos method, 
it requires less computer memory, and it works for all the 
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FIG. 2: Lanczos and Richardson results for the energies of 
the lowest zero-seniority states, relative to the non-interacting 
ground-state energy. 
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FIG. 3: Behavior of the variables near the singularity around 
the 2(^3/2 level. Variables 2:4 and xs are complex conjugates, 
xe is real over the whole range of g values. 



eigenstates, not just the ground state. Moreover, the pro- 
cedure can deal with much larger systems, well beyond 
the limits of large-scale exact diagonalizations. 

This work shows that the exactly solvable pairing mod- 
els are indeed solvable in practice, even for Fermions with 
multiple degeneracies. It opens up a whole new range 
of applications for these models. Furthermore, similar 
techniques might be useful to solve the non-linear equa- 
tions for other exactly solvable pairing and spin mod- 
els 
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